// Copyright 2020 Poofee (https://github.com/Poofee)
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
// ------------------------------------------------------------------------
/*****************************************************************************
 *                                                                           *
 *  File:    relay.h                                                         *
 *                                                                           *
 *****************************************************************************
 *                                                                           *
 *  Authors: Poofee                                                          *
 *  Email:   poofee@qq.com                                                   *
 *  Address: https://github.com/Poofee                                       *
 *  Date:    2020年07月15日                                                   *
 *                                                                           *
 *****************************************************************************/
#ifndef RELAY_H
#define RELAY_H

#include "datatype.h"
#include "slu_mt_ddefs.h"

#if !defined(ARMA_32BIT_WORD)
#define ARMA_32BIT_WORD
#endif
#include <armadillo> 

#include "meshGFace.h"
#include "meshGRegionNetgen.h"
#include "meshGRegion.h"
#include "Generator.h"
#include "GModel.h"
#include "gmsh.h"
#include "MLine.h"
#include "Field.h"

#include <vector>
#include <map>
#include <iostream>
#include <algorithm>
#include <cassert>
#include <unordered_map>

#include "tetralheral.h"

using namespace arma;

const double PI = 3.14159265358979323846;
/** 使用相对磁导率避免数字过大 **/
const double miu0 = PI * 4e-7;

template <class ForwardIterator>
ForwardIterator Myunique (ForwardIterator first, ForwardIterator last)
{
    if (first==last) return last;

    ForwardIterator result = first;
    while (++first < last)
    {
        if(*result == *first){
            /** 由于在这里已经到达last了，导致下一次加的时候就超过
             * last了，报错 **/
            while (++first != last) {
//                qDebug()<<*first;
                if(*result != *first){
                    *result = *first;
                    break;
                }
            }
        }else{
            *(++result)=*first;
//            qDebug()<<*first;
        }
        /** 最后一个数是不重复的 **/
//        if (first==last-1){
//            break;
//        }
        /** 最后几个数是重复的 **/
        if (first==last){
            return result;
        }
    }
    return ++result;
}

struct pair_hash {
	template<class T1, class T2>
	std::size_t operator() (const std::pair<T1, T2>& p) const {
		auto h1 = std::hash<T1>{}(p.first);
		auto h2 = std::hash<T2>{}(p.second);
		return h1 ^ h2;
	}
};

class Relay
{
public:
    Relay();
    ~Relay();

    virtual void init();
    /** 几何相关 **/
    void setFileName(char* fn);
	void unitConvert(double ratio);
    void openGeo();
    void moveFace(GFace *f,double dx,double dy,double dz);
    void moveRegion(GRegion *r,double dx,double dy,double dz);
    void rotateRegion(GRegion* r,double dangle,double x1, double y1,double z1,
                      double v1, double v2,double v3);
    Point3f RotateByVector(double theta, double old_x, double old_y, double old_z, 
		                                 double vx, double vy, double vz);

    /** 边界相关 **/
    void findBoundaryPoints2D(int index);
    void findBoundaryEdges2D(int index);
    void findBoundaryEdges3D(int index);
    void findBoundaryFaces3D(int index);
    /** 分网相关 **/
    void setXiantieTag(int xiantie);
    void setAirTag(int air);
	void setRemeshTag(int remesh);
    void loadMesh();
    void deleteFaceMesh(GFace* f);
    void deleteRegionMesh(GRegion* r);
    virtual void remesh(int timestep);
    void remesh(double dx, double dy);
    void remesh3DParallel(double dx, double dy, double dz);
    void remesh3DRotate(double dangle,double x1, double y1,double z1,
                                      double v1, double v2,double v3);
    void updateField();
	bool checkConnectivityMatrix();


    /** 求解器相关 **/
    void stepIncrement();
    void setPrefixName(char* fn);
    virtual void run();
	void initPhysicalVariables(int size);
    void insertSparseMatrixData(int pos, int row, int col, double values);

	void setMaterial(int index, CMaterial* mat);
	CMaterial* getMaterial(int domainIndex);

    /** 后处理相关 **/
    virtual double calcMagForce(int index, FORCEMETHOD method);
    virtual double calcMagTorque(int index, FORCEMETHOD method);
    virtual void outputResults();
    void writeVtkFile(char *fname,int num_cell,int num_type,int num_size);
	void testBoundaryResults(char fn[]);
    void printCurrentTime();

	/** 有限元分析 **/
	virtual void StaticInitial();
	virtual void TransientInitial();
	virtual void StaticAnalysis();
	virtual void TransientAnalysis();
protected:
    int num_pts;/** 读取的节点数目，需要注意的是，节点数可能包含几何节点 **/
    int num_edges;/** 棱边数 **/
    int num_ele;/** 读取的单元数目 **/
    int numDomain;/** 域数目 **/

    int num_triangle;/** 三角形单元的数目 **/
    int index_triagle;/** 三角形单元在msh文件保存的起始位置 **/

    int num_tet;/** 四面体单元的数目 **/
    int index_tet;/** 四面体单元在msh文件保存的起始位置 **/

    int tag_xiantie; /** 可移动的区域 **/
    int tag_air; /** 空气区域 **/
	int tag_remesh;/** 重分网区域，默认geo文件中该区域为最后一个区域。 **/
    int air_position;/** 在分网中，形变区域的起始位置。 **/

    int current_step;/** 当前仿真的步数 **/
    int MAX_NONLINEARSTEPS;/** 牛顿迭代法的最大迭代步数 **/
    int total_steps;/** 总共的迭代步数 **/

    int ParallelThreads;

    double  Precision;/** 计算精度 **/
    double  Relax;

    double min_time;/** 最小求解时间 **/
    double max_time;/** 最大求解时间 **/

    double min_position;/** 最小的运动位移 **/
    double max_position;/** 最大的运动位移 **/

    double unitRatio;/** 单位换算到m的比例 **/

    /** 结果 **/
    std::vector<double> timesteps;/** 每一步的时间差 **/
	std::vector<double> totaltimes;/** 时间 **/
    std::vector<double> displacements;/** 总位移 **/
    std::vector<double> Ddisplacements;/** 每一步发生的位移 **/
    std::vector<double> velocities;/** 速度 **/
    std::vector<double> accelerations;/** 加速度 **/
    std::vector<double> PhiCoil;/** 线圈磁通 **/
    std::vector<double> ICoil;/** 线圈电流 **/
    std::vector<double> magforces;/** 电磁力 **/
    std::vector<double> UCoil;/** 线圈电压 **/
    std::vector<double> IterationStep;/** 迭代步数 **/
    std::vector<double> SimulationTime;/** 仿真时间 **/

    /** 边界 **/
    std::vector<int> boundaryPoints;/** 外部边界点 **/
    std::vector<int> allPoints;/** 所有的分网点编号 **/
    std::vector<int> freePoints;/** 所有的自由分网点编号 **/
    std::vector<FEMface> boundaryFaces;/** 外部边界面 **/
    std::vector<FEMface> allFaces;/** 所有的边 **/
    std::vector<FEMedge> boundaryEdges;/** 外部边界边 **/
    std::vector<FEMedge> allEdges;/** 所有的边 **/
	std::unordered_map<std::pair<int, int>, int, pair_hash> alledgeshash;

    CNode * pmeshnode;/** 分网节点 **/
    CElement * pmeshele;/** 分网单元 **/
    CNode * pmeshnodelast;/** 上一步的分网节点 **/
    CElement * pmeshelelast;/** 上一步的分网单元 **/

    char fileName[256];/** 文件名称 **/
    char prefixName[256];/** 每一个仿真的名字 **/

    GModel* model;/** gmsh model **/
	bool RemeshOrMorph;/** true:remesh，对区域重分网;false:morph，只形变 **/

//    int		LengthUnits;/** 长度单位变量，整型 **/
    std::vector<CMaterial*> materialList;/** 定义材料类 **/
    std::map<int,CMaterial*> materialMap;/** 建立区域与材料的映射关系 **/

	umat locs;/** 为了降低二次分配时的成本 **/
	mat vals;
	vec bbJz;
	uvec node_reorder;
	uvec node_pos;
	vec bn;
	vec A;
	vec A_last;
	vec A_old;
};

#endif // RELAY_H
